% This code can be executed to determine the phenemological surface energy
% between the two tapes. Only the Peeling maxima are used. The data is 
% plotted with the x-data (xdat) being a function of the junction area and 
% the y-data (ydat) is the peak load. This will give a variable titled, 
% "fitresult," which fits the function y=mx to the data and the m is the 
% phenemological surface energy. The fit reports the value with 95% 
% confidence interval. 
% The folder titled, "Raw Data Test Files," must be open to run this code!
clear;

%Tensile Tester Trials
names=['2019-02-08_01_70.0.txt',...
    '2019-02-27_02_75.0.txt',...
    '2019-02-27_03_66.0.txt',...
    '2019-02-27_04_65.0.txt',...
    '2019-02-27_05_52.6.txt',...
    '2019-02-27_06_40.9.txt',...
    '2019-03-01_08_85.1.txt',...
    '2019-03-01_09_52.0.txt',...
    '2019-03-01_10_49.5.txt',...
    '2020-04-15_09_30.0.txt',...
    '2020-04-15_10_32.5.txt',...
    '2020-04-15_11_35.0.txt',...
    '2020-04-15_12_37.5.txt',...
    '2020-04-15_13_40.0.txt',...
    '2020-04-15_14_42.5.txt',...
    '2020-04-15_15_45.0.txt',...
    '2020-04-15_16_47.5.txt',...
    '2020-04-15_17_50.0.txt',...
    '2020-04-15_18_52.5.txt',...
    '2020-04-15_19_55.0.txt',...
    '2020-04-15_20_57.5.txt',...
    '2020-04-15_21_60.0.txt',...
    '2020-04-15_22_65.0.txt',...
    '2020-04-15_23_70.0.txt',...
    '2020-04-15_24_75.0.txt',...
    '2020-05-21_78_42.5.txt',...
    '2020-05-21_79_35.0.txt',...
    '2020-05-21_80_50.0.txt',...
    '2020-05-21_83_47.5.txt',...
    '2020-05-21_84_32.5.txt'];

for i=1:length(names)
    chunk=i+3;
    test=names(i:chunk) == '.txt';
    if all(test)
        numchar=chunk;  % number of characters in each file name
        break
    end
end
ntrials=length(names)/numchar;
jj=zeros(1,ntrials);
jj(1)=1;
kk=zeros(1,ntrials);
kk(1)=numchar;
elasticmax=[];
elasticmaxangles=[];
peelmax=[];
peelmaxangles=[];
peakload=[];
angles=[];
for i=1:ntrials
    clear Data displacement load time
    filename=names(jj(i):kk(i));
    fileID=importdata(filename);
    j=jj(i)+numchar;
    k=kk(i)+numchar;
    jj(i+1)=j;
    kk(i+1)=k;
    Data=fileID.data;
    if str2num(filename(4))==9  %older data sets require different treatment as discussed in "Figure_3B_code"
        displacement=Data(:,1).*(50E-3);  %mm
        load=abs(Data(:,2)-abs(max(Data(:,2))));  %N
        adjusttrials=[2 3 5 6 10];
        adjustdoubletroubletrials=[10];
        adjustquadtroubletrials=[];
        adjustbigtroubletrials=[2 3];
    else
        displacement=Data(:,1);  %mm
        load=Data(:,2);  %N
        adjusttrials=[13 15 17 18 19 20 22 24 81 82 85];
        adjustdoubletroubletrials=[81];
        adjustquadtroubletrials=[19];
        adjustbigtroubletrials=[];
    end
    phi=str2num(filename(15:18));
    if any(ismember(adjusttrials,str2num(filename(12:13))))
        [pks, pkloc]=findpeaks(load,displacement);
        [pks2, pkloc2]=findpeaks(pks,pkloc);
        allpks=[pks2 pkloc2];
        sortpks=sortrows(allpks,'descend');
        peelmax=[peelmax sortpks(2,1)];
        peelmaxangles=[peelmaxangles phi];
        elasticmax=[elasticmax sortpks(1,1)];
        elasticmaxangles=[elasticmaxangles phi];
        if any(ismember(adjustdoubletroubletrials,str2num(filename(12:13))))
            peelmax(length(peelmax))=sortpks(3,1);
        end
        if any(ismember(adjustquadtroubletrials,str2num(filename(12:13))))
            peelmax(length(peelmax))=sortpks(4,1);
        end
        if any(ismember(adjustbigtroubletrials,str2num(filename(12:13))))
            [pks3, pkloc3]=findpeaks(pks2,pkloc2);
            allpks=[pks3 pkloc3];
            sortpks=sortrows(allpks,'descend');
            if str2num(filename(12:13))==2
                peelmax(length(peelmax))=sortpks(2,1);
            end
            if str2num(filename(12:13))==3
                peelmax(length(peelmax))=sortpks(3,1);
            end
        end
    else
        [pks, pkloc]=findpeaks(load,displacement);
        allpks=[pks pkloc];
        sortpks=sortrows(allpks,'descend');
        peakload=[peakload sortpks(1,1)];
        angles=[angles phi];
    end
end

%Peak Load Only Trials
peakloadonlydata=importdata('Peak Load Only Raw Data.txt');
peakloadonlydata2=peakloadonlydata.data;
anglesexp1=peakloadonlydata2(1:(size(peakloadonlydata2,1)-3),1);
maxloadexp=flip(peakloadonlydata2(1:(size(peakloadonlydata2,1)-3),2))';
anglesexp=flip(anglesexp1)';


w1=19E-3;  %m
w2=12.7E-3; %m

SE1=(-2).*w1.*((cscd(anglesexp./2)).^2);
SE2=(-2).*w2.*((cscd(angles./2)).^2);
SE3=(-2).*w2.*((cscd(peelmaxangles./2)).^2);

xdat=[SE1 SE2 SE3];
ydat=[maxloadexp peakload peelmax];

hold on;

%Plot data
plot(xdat,ydat,'LineStyle','none','Marker','.','MarkerSize',15);

%Fit data to y=mx
ff1=@(m,x) m.*x;
[fitresult, gof]=fit(xdat',ydat',ff1);
fit=@(x) fitresult.m.*x;

fitresult

%Plot the fit line
xx=linspace(0.8E-4, 1.7E-4,1000);
plot(xdat,fit(xdat),'b');
plot(xx,fit(xx),'LineStyle','--','Color','b');
xlabel('-2wcsc^2(\phi/2) (m)','FontSize',14);
ylabel('Force (N)','FontSize',14);
hold off;